home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / help / filter < prev    next >
Text File  |  1994-04-25  |  3KB  |  115 lines
  1.  
  2. Synopsis:    Digital filter implementation
  3.  
  4. Syntax:        filter ( B , A , X )
  5.         filter ( B , A , X , Zi )
  6.  
  7. Description:
  8.  
  9.     Filter is an implementation of the standard difference
  10.     equation:
  11.  
  12.     y[n] = b(1)*x[n] + b(2)*x[n-1] + ... b(nb+1)*x[n-nb]
  13.              - a(2)*y[n-1] - ... a(na+1)*y[n-na]
  14.  
  15.     The filter is implemented using a method described as a
  16.     "Direct Form II Transposed" filter. More for information see
  17.     Chapter 6 of "Discrete-Time Signal Processing" by Oppenheim
  18.     and Schafer.
  19.  
  20.     The inputs to filter are:
  21.  
  22.         B    The numerator coefficients, or zeros of the
  23.             system transfer function. The coefficients are
  24.             specified in a vector like:
  25.  
  26.             [ b(1) , b(2) , ... b(nb) ]
  27.  
  28.         A    The denominator coefficients, or the poles of
  29.             the system transfer function. the coefficients
  30.             are specified in a vector like:
  31.  
  32.             [ a(1) , a(2) , ... a(na) ]
  33.  
  34.         X    A vector of the filter inputs.
  35.  
  36.         Zi    [Optional] The initial delays of the filter.
  37.  
  38.     The filter outputs are in a list with element names:
  39.  
  40.         y    The filter output. y is a vector of the same
  41.             dimension as X.
  42.  
  43.         zf    A vector of the final values of the filter
  44.             delays.
  45.  
  46.     The A(1) coefficient must be non-zero, as the other
  47.     coefficients are divided by A(1).
  48.  
  49.     Below is an implementation of filter() in a r-file - it is
  50.     provided for informational usage only.
  51.  
  52.     #--------------------------------------------------------------
  53.     #  Simplistic version of RLaB's builtin function filter()
  54.     #  Y = filter ( b, a, x )
  55.     #  Y = filter ( b, a, x, zi )
  56.     #
  57.     
  58.     rfilter = function ( _b , _a , x , zi )
  59.     {
  60.       local (M, N, NN, a, b, k, n, ntotal, v, vi, y)
  61.     
  62.       # Copy _a and _b cause we are going to modify them
  63.       a = _a;
  64.       b = _b;
  65.     
  66.       ntotal = x.nr * x.nc;
  67.       M = b.nr * b.nc;
  68.       N = a.nr * a.nc;
  69.       NN = max ([ M, N ]);
  70.       y = zeros (x.nr, x.nc); 
  71.     
  72.       # Fix up pole and zero vectors.
  73.       # Make them the same length, this makes
  74.       # filter's job much easier.
  75.     
  76.       if (N < NN) { a[NN] = 0; }
  77.       if (M < NN) { b[NN] = 0; }
  78.       
  79.       # Adjust filter coefficients
  80.       if (a[1] == 0) { error ("rfilter: 1st A term must be non-zero"); }
  81.       a[2:NN] = a[2:NN] ./ a[1];
  82.       b = b ./ a[1];
  83.     
  84.       # Create delay vectors and load inital delays.
  85.       # Add an extra term to vi[] to make filter's 
  86.       # job a little easier. This extra term will
  87.       # always be zero.
  88.     
  89.       v = zeros (NN, 1);
  90.       vi = zeros (NN+1, 1);
  91.     
  92.       if (exist (zi))
  93.       {
  94.         vi[1:NN] = zi;   
  95.       }
  96.     
  97.       #
  98.       # Do the work...
  99.       #
  100.     
  101.       for (n in 1:ntotal)
  102.       {
  103.         v[1] = b[1]*x[n] + vi[2];
  104.         y[n] = v[1];
  105.         for (k in 2:NN)
  106.         {
  107.           v[k] = b[k]*x[n] - a[k]*v[1] + vi[k+1];
  108.           vi[k] = v[k];
  109.         }
  110.       }
  111.     
  112.       return << y = y; zf = v >>;
  113.     };
  114.     #--------------------------------------------------------------
  115.